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ABSTRACT 

We perform kinetic simulations of diffusive shock acceleration (DSA) in Type la supernova remnants (SNRs) 
expanding into a uniform interstellar medium (ISM). Bohm-like diffusion due to self-excited Alfven waves is as- 
sumed, and simple models for Alfvenic drift and dissipation are adopted. Phenomenological models for thermal 
leakage injection are considered as well. We find that the preshock gas temperature is the primary parameter that 
governs the cosmic ray (CR) acceleration efficiency and energy spectrum, while the CR injection rate is a secondary 
parameter. For SNRs in the warm ISM of To <, lO^K, if the injection fraction is ^ ^ 10-"*, the DSA is efficient 
enough to convert more than 20 % of the SN explosion energy into CRs and the accelerated CR spectrum exhibits 
a concave curvature flattening to E~^-^, which is characteristic of CR modified shocks. Such a fiat source spectrum 
near the knee energy, however, may not be reconciled with the CR spectrum observed at Earth. On the other 
hand, SNRs in the hot ISM of Tq « lO^K with a small injection fraction, ^ < 10"'*, are inefficient accelerators with 
less than 10 % of the explosion energy getting converted to CRs. Also the shock structure is almost test-particle 
like and the ensuing CR spectrum can be steeper than E~'^. With amplified magnetic field strength of order of 
30/iG, Alfven waves generated by the streaming instability may drift upstream fast enough to make the modified 
test-particle power-law as steep as E^^'^, which is more consistent with the observed CR spectrum. 

Key Words : cosmic ray acceleration - supernova remnants ~ hydrodynamics - methods:numerical 



I. INTRODUCTION 

It is believed that most of the Galactic cosmic rays 
(CRs) are accelerated in the blast waves driven by 
supernova (SN) explosions (e.g., Blandford & Eichler 
1987, Reynolds 2008 and references therein). If about 
10 % of Galactic SN luminosity, Lsn ~ lO'^^erg s~*, 
is transfered to the CR component, the diffusive shock 
acceleration (DSA) at supernova remnants (SNRs) can 
provide the CR luminosity, Lcr ~ lO^'^erg s~^ that es- 
capes from the Galaxy. Several time-dependent, kinetic 
simulations of the CR acceleration at SNRs have shown 
that an order of 10 % of the SN explosion energy can be 
converted to CRs, when a fraction ~ 10~^ of incoming 
thermal particles are injected into the CR population 
at the subshock {e.g., Berezhko, & Volk 1997; Berezhko 
et al 2003; Kang 2006). 

X-ray observations of young SNRs such as SN1006 
and RCW86 indicate the presence of 10-100 TeV elec- 
trons emitting nonthermal synchrotron emission imme- 
diately inside the outer SNR shock (Koyama et al. 1995; 
Bamba et al. 2006, Helder et al. 2009). They provide 
clear evidence for the efficient acceleration of the CR 
electrons at SNR shocks. Moreover, HESS gamma- 
ray telescope detected TeV emission from several SNRs 
such as RXJ1713. 7-3946, Cas A, Vela Junior, and 
RCW86, which may indicate possible detection of 7r° 
7- rays produced by nuclear collisions of hadronic CRs 
with the surrounding gas (Aharonian et al. 2004, 2009; 
Berezhko & Volk 2006; Berezhko et al. 2009; Morhno 
et al. 2009, Abdo et al. 2010). It is still challenging to 
discern whether such emission could provide direct evi- 



dence for the acceleration of hadronic CRs, since 7-ray 
emission could be produced by inverse Compton scat- 
tering of the background radiation by X-ray emitting 
relativistic electrons. More recently, however, Fermi 
LAT has observed in GeV range several SNRs inter- 
acting with molecular clouds, providing some very con- 
vincing evidence of 7r° decay 7-rays (Abdo et al. 2009, 
2010). 

In DSA theory, a small fraction of incoming thermal 
particles can be injected into the CR population, and 
accelerated to very high energies through their inter- 
actions with resonantly scattering Alfven waves in the 
converging flows across the SN shock {e.g., Drury et 
al. 2001). Hence the strength of the turbulent mag- 
netic field is one of the most important ingredients, 
which govern the acceleration rate and in turn the max- 
imum energy of the accelerated particles. If the mag- 
netic field strength upstream of SNRs is similar to the 
mean interstellar medium (ISM) field of i?isM ^ 5/iG, 
the maximum energy of CR ions of charge Z is esti- 
mated to be -Emax ^ eV (Lagage & Cesarsky 
1983). However, high- resolution X-ray observations of 
several young SNRs exhibit very thin rims, indicating 
the presence of magnetic fields as strong as a few 100/iG 
downstream of the shock {e.g., Bamba et al. 2003, Pari- 
zot et al. 2006). Moreover, theoretical studies have 
shown that efficient magnetic field amplification via 
resonant and non-resonant wave-particle interactions 
is an integral part of DSA (Lucek & BeU 2000, Bell 
2004). If there exist such amplified magnetic fields in 
the upstream region of SNRs, CR ions might gain en- 
ergies up to i?max ^ lO^^'^Z eV, which may explain 
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the all-particle CR spectrum up to the second knee at 
~ 10^^ cV with rigidity-dependent energy cutoffs. A 
self-consistent treatment of the magnetic field amplifi- 
cation has been implemented in several previous stud- 
ies of nonlinear DSA {e.g., Amato & and Blasi 2006, 
Vladimirov et al. 2008). 

In Kang 2006 (Paper I, hereafter), we calculated 
the CR acceleration at typical remnants from Type 
la supernovae expanding into a uniform interstellar 
medium (ISM). With the upstream magnetic fields of 
Bq = 30/LiG amplified by the CR streaming instabil- 
ity, it was shown that the particle energy can reach 
up to lO^^Z eV at young SNRs of several thousand 
years old, which is much higher than what Lagage & 
Cesarsky predicted. But the CR injection and accel- 
eration efficiencies are reduced somewhat due to faster 
Alfven wave speed. With the particle injection frac- 
tion ~ 10"^ - 10-3, the DSA at SNRs is very efficient, 
so that up to 40-50 % of the explosion energy can 
be transferred to the CR component. We also found 
that, for the SNRs in the warm ISM (To = lO^'K), 
the accelerated CR energy spectrum should exhibit a 
concave curvature with the power-law slope, a (where 
N{E) cx E-°') flattening from 2 to 1.6 at E > 0.1 
TeV. In fact, the concavity in the CR energy spec- 
trum is characteristic of strong (M > 10) CR modi- 
fied shocks when the injection fraction is greater than 
10"^. {e.g., Malkov & Drury 2001, Berezhko & Volk 
1997, Blasi et al. 2005) 

Recently, Ave et al. (2009) have analyzed the spec- 
trum of CR nuclei up to ~ 10^'' eV measured by 
TRACER instrument and found that the CR spec- 
tra at Earth can be fitted by a single power law of 
J{E) (X E~'^'^'^ . Assuming an energy-dependent prop- 
agation path length (A cx i?^°-^), they suggested that 
a soft source spectrum, N{E) with a ~ 2.3 — 2.4 is pre- 
ferred by the observed data. However, the DSA pre- 
dicts that a = 2.0 for strong shocks in the test-particle 
limit and even smaller values for CR modified shocks 
in the efficient acceleration regime as shown in Paper 
I. Thus in order to reconcile the DSA prediction with 
the TRACER data the CR acceleration efficiency at 
typical SNRs should be minimal and perhaps no more 
than 10 % of the explosion energy transferred to CRs 
{i.e., test-particle limit). Moreover, recent Fermi-LAT 
observations of Cas A, which is only 330 years old and 
has just entered the Sedov phase, indicate that only 
about 2% of the explosion energy has been transfered 
to CR electrons and protons, and that the soft proton 
spectrum with E~^'^ is preferred to fit the observed 
gamma-ray spectrum (Abdo et al. 2010). According 
to Paper I, such inefficient acceleration is possible only 
for SNRs in the hot phase of the ISM and for the in- 
jected particle fraction smaller than 10-**. One way 
to soften the CR spectrum beyond the canonical test- 
particle slope {a > 2) is to include the Alfvenic drift in 
the precursor, which reduces the velocity jump across 
the shock. Zirakashvili & Ptuskin (2008) showed that 
the Alfvenic drift in the amplified magnetic fields both 



upstream and downstream can drastically soften the 
accelerated particle spectrum. We will explore this is- 
sue using our numerical simulations below. 

Caprioli et al. (2009) took a different approach to 
reconcile the concave CR spectrum predicted by non- 
linear DSA theory with the softer spectrum inferred 
from observed J{E). They suggested that the CR spec- 
trum at Earth is the sum of the time integrated flux of 
the particles that escape from upstream during the ST 
stage and the flux of particles confined in the remnant 
and escaping at later times. They considered several 
cases and found the injected spectrum could be softer 
than the concave instantaneous spectrum at the shock. 
The main uncertainties in their calculations are related 
with specific recipes for the particle escape. It is not 
well understood at the present time how the particles 
escape through a free escape boundary (xesc) located 
at a certain distance upstream of the shock or through 
a maximum momentum boimdary due to lack of (self- 
generated) resonant scatterings above an escape mo- 
mentum. The escape or release of CRs accelerated in 
SNRs to the ISM remains largely unknown and needs 
to be investigated further. 

One of the key aspects of the DSA model is the 
injection process through which suprathermal parti- 
cles in the Maxwellian tail get accelerated and injected 
into the Fermi process. However, the CR injection 
and consequently the acceleration efficiency still re- 
main uncertain, because complex interplay among CRs, 
waves, and the underlying gas flow {i.e., self-excitation 
of waves, resonant scatterings of particles by waves, 
and non-linear feedback to the gas flow) is all model- 
dependent and not understood completely. 

In this paper, we adopted two different injection 
recipes based on thermal leakage process, which were 
considered previously by us and others. Then we have 
explored the CR acceleration at SNR shocks in the dif- 
ferent temperature phases {i.e., different shock Mach 
numbers) and with different injection rates. Details of 
the numerical simulations and model parameters are 
described in §11. The simulation results are presented 
and discussed in §111, followed by a summary in §IV. 

II. NUMERICAL METHOD 
(a) Spherical CRASH code 

Here we consider the CR acceleration at a quasi- 
parallel shock where the magnetic fleld lines are paral- 
lel to the shock normal. So we solve the standard gas- 
dynamic equations with CR pressure terms added in 
the Eulerian formulation for one dimensional spherical 
symmetric geometry. The basic gasdynamic equations 
and details of the spherical CRASH (Cosmic- Ray Amr 
SHock) code can be found in Paper I and Kang & Jones 
(2006). 

In the kinetic equation approach to numerical study 
of DSA, the following diffusion-convection equation for 
the particle momentum distribution, f{p), is solved 
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along with suitably modified gasdynamic equations 
{e.g., Kang & Jones 2006): 



^2 



2 / N^S 

r K{r,y) — 



(1) 



where g = j//, with f{p,r,t) the pitch angle averaged 
CR distribution, and y = ln(p), and K(r, y) is the diffu- 
sion coefficient parallel to the field lines . So the proton 
number density is given by ncR.p = 4:11 J f(p)p'^dp. For 
simplicity we express the particle momentum, p in units 
of rUpC and consider only the proton component. 

The velocity represents the effective relative mo- 
tion of scattering centers with respect to the bulk flow 
velocity, u. The mean wave speed is set to the Alfvcn 
speed, i.e., Uw = va = B/^Airp in the upstream re- 
gion. This term reflects the fact that the scattering 
by Alfven waves tends to isotropize the CR distribu- 
tion in the wave frame rather than the gas frame. In 
the postshoek region, = is assumed, since the 
Alfvenic turbulence in that region is probably relatively 
balanced. This reduces the velocity difference between 
upstream and downstream scattering centers compared 
to the bulk flow, leading to less efficient DSA. This in 
turn affects the CR spectrum, and so the 'modifled' 
test-particle slope can be estimated as 



9tp = 



3(mo - va) 

Uo- VA - U2 



(2) 



where f{p) oc p^"?'? is assumed {e.g., Kang et al. 2009). 
Hereafter we use the subscripts '0', '1', and '2' to de- 
note conditions far upstream of the shock, immedi- 
ately upstream of the gas subshock and immediately 
downstream of the subshock, respectively. Thus the 
drift of Alfven waves in the upstream region tends to 
soften the CR spectrum from the canonical test-particle 
spectrum of /(p) oc p~^ if the Alfven Mach number 
{Ma = Us/va) is small. We note a = g — 2 for rela- 
tivistic energies. For example, for a strong shock with 
U2 = uo/4 in the test particle limit, we can obtain the 
observed value of a = 2.3 if = 0.173uo. 

Gas heating due to Alfven wave dissipation in the 
upstream region is represented by the term 



W{r,t) = -ujhVa 



dPc 
dr ' 



(3) 



where Pc = {AnmpC^/S) J g{p)dp/yjp^ + 1 is the CR 
pressure. This term is derived from a simple model 
in which Alfven waves are amplified by streaming CRs 
and dissipated locally as heat in the precursor region 
{e.g., Jones 1993). As was previously shown in SNR 
simulations {e.g., Berezhko & Volk 1997, Kang & Jones 
2006), precursor heating by wave dissipation reduces 
the subshock Mach number thereby reducing DSA ef- 
ficiency. The parameter wjj is introduced to control 



the degree of wave dissipation. We set ujh = 1 for all 
models unless stated otherwise. 

Accurate solutions to the CR diffusion-convection 
equation require a computational grid spacing signifi- 
cantly smaller than the particle diffusion length. Ace <C 
Xd{p) = k{p)/us. With Bohm-like diffusion coefficient, 
K,{p) OC p, a wide range of length scales must be re- 
solved in order to follow the CR acceleration from the 
injection energy (typically piaj ^ 10~^) to highly rela- 
tivistic energy {p ^ 1). This constitutes an extremely 
challenging numerical task, requiring rather extensive 
computational resources. In order to overcome this dif- 
ficulty, we have developed CRASH code in ID plane- 
parallel geometry (Kang et al. 2001) and in ID spheri- 
cal symmetric geometry (Kang & Jones 2006) by com- 
bining Adaptive Mesh Refinement technique and sub- 
grid shock tracking technique. Moreover, we solve the 
fluid and diffusion-convection equations in a frame co- 
moving with the outer spherical shock in order to im- 
plement the shock tracking technique effectively in an 
expanding spherical geometry. In the comoving grid, 
the shock remains at the same location, so the compres- 
sion rate is applied consistently to the CR distribution 
at the subshock, resulting in much more accurate and 
efficient low energy CR acceleration. 

(b) Injection Recipes for Thermal Leakage 

The injection rate with which suprathermal particles 
are injected into CRs at the subshock depends in gen- 
eral upon the shock Mach number, field obliquity angle, 
and strength of Alfven turbulence responsible for scat- 
tering. In thermal leakage injection models suprather- 
mal particles well into the exponential tail of the post- 
shock Maxwellian distribution leak upstream across a 
quasi-parallel shock (Malkov & Volk 1998; Malkov & 
Drury 2001). Currently, however, these microphysics 
issues are known poorly and any quantitative predic- 
tions of macrophysical injection rate require extensive 
understandings of complex plasma interactions. Thus 
this process has been handled numerically by adopting 
some phenomenological injection schemes in which the 
particles above a certain injection momentum pinj cross 
the shock and get injected to the CR population. 

There exist two types of such injection models con- 
sidered previously by several authors. In a simpler 
form, pinj represents the momentum boundary between 
thermal and CR population and so the particles are in- 
jected at this momentum {e.g., Kang & Jones 1995, 
Berezhko & Volk 1997, Blasi et al. 2005). The injec- 
tion momentum is then expressed as 



Pinj = RinjPth, 



(4) 



where R-mj is a constant and pth = \/'2kBT2mp is the 
thermal peak momentum of the Maxwellian distribu- 
tion of the immediate postshoek gas with temperature 
T2, and ks is the Boltzmann constant. The CR distri- 
bution at pinj is then fixed by the Maxwellian distribu- 
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Table 1. Model Parameters 



ModcP 


UH (ISM) 


To 


Eo 




ro 


to 


Uo 


Po 




(cm""*) 


(K) 


(10^^ ergs) 


(/iG) 


(pc) 


(years) 


(10^ km s-i) 


(10-^erg cm-3) 


WA/WB 


0.3 


3.3 X 10'' 


1. 


30 


3.19 


255. 


1.22 


1.05 


MA/MB 


0.03 


10^ 


1. 


30 


6.87 


549. 


1.22 


1.05 X 10-1 


HA/HB 


0.003 


10^ 


1. 


30 


14.8 


1182. 


1.22 


1.05 X 10-2 



^ 'W, 'M', and 'H' stands for the warm, intermediate, and hot phase of the ISM, respectively, while 'A' and 'B' 
stands for the injection recipes A and B, respectively, described in II (b). 



tion, 

/(Pinj) =n2 i 27^^3X2 ) ^"^^^'^^i^' 

where n2 is the postshock proton number density. Thus 
the constant parameter i?inj controls the injection rate 
in this model. Here wc refer this as 'injection recipe B' 
and consider the cases of i?inj = 3.6 and 3.8. 

In Kang et al. (2002), on the other hand, a smooth 
"transparency function", rosc(eB,w) is adopted, rather 
than a step-like filter function of the injection recipe 
B. This function expresses the probability of supra- 
thermal particles at a given velocity, v, leaking up- 
stream through the postshock MHD waves. One free 
parameter controls this function; es = Bo/B±, which 
is the inverse ratio of the amplitude of the postshock 
MHD wave turbulence B± to the general magnetic field 
aligned with the shock normal, Bq (Malkov & Volk 
1998). In this model, the leakage probability Tgsc > 
above pi ~ mpU2(l + l.OT/es) oc pth, and the "effec- 
tive" injection momentum is a few times pi. So the 
injection momentum can be expressed as 

Pinj = Qinj {Ms , eB)Pth ■ (6) 

Note that the ratio Qinj is a function of the subshock 
Mach number, Mg, as well as the parameter eb, while 
the constant ratio i?inj is independent of Ms- The value 
of Qinj is larger (and so the injection rate is smaller) 
for weaker subshocks and for smaller cb (see Kang et 
al. 2002). In an evolving CR modified shock, the sub- 
shock weakens as the precursor develops due to non- 
linear feedback of the CR pressure and so the injection 
rate decreases in time. We refer this as 'injection recipe 
A' and consider 0.2 < < 0.3 here. 

In Paper I we only considered the gas with pro- 
tons (i.e., mean molecular weight fi = 1), but here 
we assume fully ionized plasma with cosmic abundance 
{fi = 0.61). As a result, for given gas pressure and den- 
sity, the temperature is lower and so slightly larger cb 
is needed to obtain the similar level of injection as in 
Paper I. Note that es = 0.16 - 0.2 in Paper I. 

The efficiency of the particle injection is quantified 
by the fraction of particles swept through the shock 



that have been injected into the CR distribution: 

_ /47rr2dr/47r/(p,r,t)p2dp 
' /47rr2noU,dt ' ^ ' 

where no is the proton number density far upstream 
and Ts is the shock radius. Recent observations of non- 
thermal radiation from several SNRs indicate that the 
injection fraction is about ^ ~ IQ-"' {e.g., Berezhko et 
al. 2009, Morhno et al. 2009). 

In our simulations, initially there is no pre-existing 
CRs and so all CR particles are freshly injected at the 
shock. 

(c) A Bohm-Iike Diffusion Model 

Self-excitation of Alfven waves by the CR streaming 
instability in the upstream region is an integral part of 
the DSA (Beh 1978; Lucek & Bell 2004). The particles 
are resonantly scattered by those waves, diffuse across 
the shock, and get injected into the Fermi first-order 
process. These complex interactions are represented 
by the diffusion coefficient, which is expressed in terms 
of a mean scattering length. A, and the particle speed, 
V, as k{x,p) = Xv/3. The Bohm diffusion model is 
commonly used to represent a saturated wave spectrum 
{i.e., \ = rg, where rg is the gyro-radius), kb{p) = 
i^nP^I {f + 1)^/^- Here /t„ = mc^/{2,eB) = 3.13 x 
lO^^cm^s-ifi-i, and B^^ is the magnetic field strength 
in units of microgauss. As in Paper I, we adopt a Bohm- 
like diffusion coefficient that includes a weaker non- 
relativistic momentum dependence, 

K{r,p) = Kn-P-^. (8) 

Since we do not follow explicitly the amplification of 
magnetic fields due to streaming CRs, we simply as- 
sume that the field strength scales with compression 
and so the diffusion coefficient scales inversely with den- 
sity. 

III. Simulations of Sedov- Taylor Blast Waves 
(a) SNR Model Parameters 

As in Paper I, we consider a Type la supernova ex- 
plosion with the ejecta mass, Mej = 1.4M0, expand- 
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Fig. 1. — Time evolution of SNR model MA with es = 0.25 (upper panels) and SNR model MB with Rinj = 3.6 (lower 
panels) at t/to = 1., 3., 6., 10. and 15. In the right panels, heavy lines are for the CR pressure, while thin lines are for the 
gas pressure. The model parameters are Mej — 1.4Mq, Eo = 10^^ ergs, nn = 0.03cm~'^, To — lO^K, and Bo — 30^G. See 
Table 1 for the normalization constants. 



ing into a uniform ISM. All models have the explo- 
sion energy, Eq — 10^^ ergs. Previous studies have 
shown that the shock Mach number is the key param- 
eter determining the evolution and the DSA efheiency, 
although other processes such as the particle injec- 
tion, the Alfvenic drift and dissipation do play certain 
roles (e.g. Kang & Jones 2002, 2007). So here three 
phases of the ISM are considered: the warm phase with 
uh — 0.3cm~^ and Tq = 3 x lO^K, the hot phase with 
uh = 0.003 cm""^ and Tq — lO^K, and the intermediate 
phase with nn = 0.03 cm"^ and To = lO^K. The pres- 
sure of the background gas is Pism ~ 10~^^ erg cm~^. 
Model parameters are summarized in Table 1. For ex- 
ample, 'WA' model stands for the warm phase and the 
injection recipe A, while 'MB' model stands for the in- 
termediate phase and the injection recipe B. 

Recent X-ray observations of young SNRs indicate 
a magnetic field strength much greater than the mean 



ISM field of 5^G {e.g., Berezhko et al. 2003; Volk et 
al. 2005). Thus, to represent this effect we take the 
upstream field strength, Bq = 30/xG. The strength of 
magnetic field determines the size of diffusion coefR- 
cient, K„, and the drift speed of Alfven waves relative 
to the bulk flow. The Alfven speed is given by va = 
VA,id{p/ Po)~^/'^ where wa,o = (1-8 kms"i)B^,/^ni". So 
in the hot phase of the ISM (HA/HB models), va^ = 
986 kms^^ and va,q/us ~ 0.175 a.t t = to, 

The physical quantities are normalized, both in the 
numerical code and in the plots below, by the following 
constants: 



3M, 



to 



Airpo 
Eo 



1/3 



1/2 
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WA: 



MA: ejj=0.25 



HA: ejj=0.25 




1 10 1 10 1 10 

t/t„ t/t„ t/t. 

Fig. 2. — Immediate pre-subshock density, pi, post-subshock density, p2, post-subshock CR and gas pressure in units of 
the ram pressure of the unmodified Sedov- Taylor solution, poUgj^ oc (t/to)~''^^ , the CR injection parameter, ^, and subshock 
Mach number, Mg are plotted for models WA (left panels), MA (middle panels), and HA (right panels). See Table 1 for the 
model parameters. The injection recipe A with es = 0.25 is adopted. For WA and HA models, the dashed lines are for the 
runs with a reduced wave heating parameter, wh = 0.5. 



Po = (2.34 X 10-^^gcm-^)nH, 

Po = Poul- 

These values are also given in Table 1 for reference. 
Note that these physical scales depend only on uh, 
since M^j and Eq are the same for all models. 

For a SNR propagating into a uniform ISM, the high- 
est momentum, pmax, is achieved at the beginning of 
the Sedov- Taylor (ST hereafter) stage and the trans- 
fer of explosion energy to the CR component occurs 
mostly during the early ST stage {e.g., Berezhko et 
al. 1997). In order to take account of the CR acceler- 
ation from free expansion stage through ST stage, we 
begin the calculations with the ST similarity solution 
at t/to = 0.5 and terminated them at t/to = 15. See 
Paper I for further discussion on this issue. 



IV. RESULTS 

(a) Remnant Evolution 

Fig. 1 shows the evolution of SNRs in the intermedi- 
ate temperature phase with es = 0.23 (injection recipe 
A) and with i?inj — 3.6 (injection recipe B). The spa- 
tial profile and the evolution of SNRs are quite similar 
in the two models, implying that detail difference be- 
tween the two injection recipes is not crucial. This is 
because the subshock Mach number reduces to Ms « 4 
at the early stage and remains the same until the end of 
simulations in the both models, resulting in similar evo- 
lutionary behavior of f (see Fig. 2). In these efficient 
acceleration models, by the early ST stage, t/to ~ 1, 
the forward shock has already become dominated by 
the CR pressure and the total density compression ra- 
tio becomes P2//O0 ~ 6 in both models. Spatial distri- 
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WA: £^ = 0.2 3 



MA: £b = 0.2 0.3 



HA: ejj=0.2, 0.3 




1 10 1 10 1 10 

t/t„ t/t„ t/t. 

Fig. 3. — The same as Fig. 2 except that es = 0.2 (solid lines) or 0.3 (dashed lines) is adopted 



bution of the CR pressure widens and becomes broader 
than that of the gas pressure at the later stage. 

The precursor length scale is given by the diffusion 
length of the highest momentum particles, /d.max ~ 
0.1 J Us{t)dt, independent of the diffusion coefficient k„ 
(Kang et al. 2009). In the test-particle limit, the shock 
would follow the ST similarity solution given by 



Ust/uo = OAUt/toY 



(9) 



where = 1.15167 is the similarity parameter (Spritzer 
1978). Then /d,max/ro ~ O.l^,{t/to)"-'^. Since the 
shock radius of the ST solution is tst/^'o = ^s{t/to)^''^, 
the relative width of the precursor is estimated to be 
^d,max/''s ~ 0.1, independent of k„. In Fig. 1 one can 
see narrow precursors in the density and CR pressure 
distribution, consistent with this estimate. 

The shock Mach number is the primary parameter 
that determines the CR acceleration efficiency, while 
the injection parameters es and i?inj are the secondary 
parameters. So the temperature (i.e., sound speed) of 



the background ISM is important. The Mach numbers 
of the initial shock are M^^j « 310, 180, and 60 in the 
warm, intermediate, hot ISM models, respectively. For 
the warm (WA/WB) and intermediate (MA/MB) mod- 
els, the CR acceleration is efficient with the postshock 
CR pressure, 0.2 < Pc,2/(poC^|t) ^ 0-4, and the shock 
is significantly modified. We will refer these models as 
'efficient acceleration models'. For HA/HB models, the 
CR acceleration is inefficient with -Pc,2/(poC^|t) < 0-1 
and the shock is almost test-particle like. So the hot 
ISM models are 'inefficient acceleration models'. Re- 
garding the injected particle fractions, the models with 
es > 0.23 or i?i„j — 3.6 — 3.8 have the injection frac- 
tion, ^ > 10~^ and represent 'efficient injection mod- 
els'. The models with es — 0.2 have ^ « lO^**-^ and 
are 'inefficient injection models'. 

Figs. 2-4 show the evolution of shock properties such 
as the compression factors, postshock pressures, the in- 
jection fraction, and subshock Mach number for vari- 
ous models. In Fig. 2, the models with — 0.25 are 
shown for ljh = 0.5 (dashed lines) and ljh = 1-0 (solid 
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Fig. 4. — The same as Fig. 2 except that the injection recipe B with i?inj = 3.6 or iiinj — 3.8 is adopted 



lines). We can see that the more efficient wave dissipa- 
tion (i.e., larger ujh) reduces the CR acceleration and 
the flow compression. Here WA and MA models are 
characterized by both efficient injection and efficient 
acceleration, while HA models show an efficient injec- 
tion but inefficient acceleration. Most of shock proper- 
ties seem to approach to more or less time-asymptotic 
values before t/to = 1. As the precursor grows, the 
subshock weakens to 3 < Ms < 5 in these models. The 
injected CR particle fraction is about ^ « 10^"^. The 
postshock CR pressure is about Pc,2 / {poUgrp) w 0.4, 
0.25, and 0.1 for WA, MA, and HA models, respec- 
tively. The compression factor in the precursor varies 
a little among different models, typically pi/po ~ 2 — 3. 
The total density compression is P2/P0 ~ 7-8 for WA 
models, 5.5 for MB model, 4.5 for HA models, indicat- 
ing the CR modified shock structure. 

Comparison of Figs. 2 and 3 tells us how the 
CR acceleration depends on the injection parameter 
cb and consequently on ^. The injection fraction is 
^ w 10~^'^ for €b = 0.2 (inefficient injection mod- 



els), C w 10-3 - 10-3-5 for 0.23 < es < 0.25, and 
^ « IQ-^-^ — 10^^ for eB = 0.3. In the inefficient injec- 
tion models, the subshock Mach number and postshock 
properties change rather gradually throughout the Se- 
dov stage and the shock is almost test-particle like with 
Pi/po ~ 1- On the other hand, for = 0.3 the injec- 
tion fraction seems much higher than what is inferred 
from recent observations of nonthermal emission from 
young SNRs {e.g., Morhno et al. 2009). The postshock 
CR pressure, Pc,2 / [poUgj,) , is roughly independent of 
the injection fraction as long as ^ > lO"^-^ (e^ ^ 0.23). 
But, for inefficient injection models with — 0.2, this 
ratio is reduced significantly. 

We can see that the models with = 0.25 (Fig. 2) 
and the models with = 3.6 (Fig. 4) have similar 
results. This confirms that the shock Mach number is 
the primary factor that controls the CR acceleration. 
For the injection recipe B, the injection rate does not 
depend on the subshock Mach number, so the evolu- 
tion of ^ is similar among WB, MB, and HB models. 
The models with _Rinj = 3.8 have about 3 times smaller 
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Fig. 5. — The CR distribution at tlie slioclc, g(rs,p), and its slope, q{p) = —d{lng{rs,p))/d\np + 4, tlie volume integrated 
CR number, G{p) = J g{r,p)ATvr^dr, and its slope, Q{p) = — d(lnG(p))/cilnp -I- 4, are shown at t/to — 1, 3, 6, 10, and 15 
for models WA (left panels), MA (middle panels), and HA (right panels). The injection recipe A with en = 0.25 is adopted. 



injection fraction and evolve more slowly, compared to 
those with = 3.6. But both models seems to ap- 
proach to the similar states at t/to > 10. 

(b) Cosmic Ray Spectrum 

The rate of momentum gain for a particle under go- 
ing DSA is given by 

-n = — o — { — + — ]p- y^^) 

at 6 \Kq K2 J 

Assuming that the shock approximately follows the ST 
solution and that the compression ratio is about four, 
the maximum momentum of protons accelerated from 
ti to t, can be calculated as 



(11) 



For our simulations started from ti/to = 0.5, this 
asymptotes to Pmax ~ O-^^iUgto/ Kn) ~ 10^-^ at large t, 
which corresponds to ii^max ~ 10^^-^ eV. 

Figs. 5-7 show the CR distribution function at the 
shock, g{rs,p), and its slope q{p) = —d{lTig{rs,p))/dliip+ 
4, and the volume integrated CR spectrum, G{p) = 
J 4:TTg{r,p)r'^dr audits slope (5(q) = —d{\nG{p))/d\n.p+ 
4. The thermal population is included in the plot of 
g{rs,p) in order to demonstrate how it is smoothly 
connected with the CR component through thermal 
leakage. For the volume integrated spectrum, only the 
CR component is shown. We note that in our sim- 
ulations the particles escape from the simulation box 
only by diffusing upstream and no escape condition is 
enforced. Thus G{p) represents the spectrum of the 
particles confined by the shock, including the particles 
in the upstream region. From the spectra in Figs. 5-7, 
we can see that Pmax approaches up to ^ 10^^ — 10^^-^ 
eV/c at t/to = 15 for all models, which is consistent 
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with the estimate given in Eq. (11). 

In Fig. 5, the CR spectra at the shock in the models 
with es = 0.25 (high injection rate with ^ « 10~^) ex- 
hibit the canonical nonlinear concave curvature. This 
is a consequence of the following two effects: the large 
compression factor across the shock structure and the 
decreasing injection rate due to the slowing of the shock 
speed. With the CR modified flow structure, the slope 
near Pmax becomes harder with qt = 3sat/{sat — 1), 
where s = 1 — va/u^ is the modification factor due 
to the Alfvenic drift and at — P2/P0 ^ 4 is the total 
shock compression ratio. If the subshock Mach num- 
ber reduces to Ms « 3 — 5, the test-particle slope at 
low momenta becomes Qs = 2>sa s / {sa s — 1) w 4.2 — 4.5, 
where as — P2I P\ is the compression ratio across the 
subshock. The particle flux through the shock, poUs, 
decreases, because the SNR shock slows down in time. 
At the same time the injection rate decreases because 
the injection process is less efficient at weaker shocks. 
The combined effects result in the reduction of the am- 
plitude of f{rs,p) near pi^y Consequently, the CR 



spectrum at lower momentum steepens and decreases 
in time. Fig. 5 demonstrates that the modified flow 
structure along with slowing down of the shock speed 
accentuates the concavity of the CR spectrum in much 
higher degrees than what is normally observed in plane- 
parallel shocks. 

However, the volume integrated spectrum G(p) is 
more relevant for unresolved observations of SNRs or 
for the total CR spectrum produced by SNRs. The 
concavity of G{p) is much less pronounced than that 
of g{rs,p), and its slope Q{p) varies 4.2-4.4 at low mo- 
mentum and 3.5-4.0 at high momentum among differ- 
ent models. We see that G{p) stays almost constant for 
t/to > 6, especially for 10" < E < lO^^eV, while ex- 
tending in the momentum space with decreasing Pmin 
and increasing Pmax- This can be understood as fol- 
lows. From Figs. 2-4, Pc.2/(poC^st) ^constant for 
t/to > 1 (except in the inefficient injection model with 
e_B — 0.2), so the CR pressure evolves like Pc,2 oc t~^^^ 
(see also Fig. 1). The total CR energy associated 
with the remnant is roughly -EcR oc Pc^st constant. 
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Fig. 7. — The same as Fig. 5 except that the injection recipe B with i?inj = 3.8 is adopted. 



Since Eck cx: J G{p)dp approximately, so G{p) should 
approach to an asymptotic form for t/to ^ 1. In other 
words, the distribution function g{r/rs,p) decreases as 
^-6/5 ^gj-mg of the normalized coordinate, r/r^, but 
the volume occupied by the remnant increases as t^^^, 
resulting in more or less constant G{p). Using this and 
the fact that Pmax asymptotes to 0.6{ulto/ Kn) at large 
t, we can predict that the form of G{p) would remain 
about the same at much later time. 

As discussed in the Introduction, the CR spectrum 
observed at Earth is J{E) cx E''^ '^'' for 10^ < E < lO" 
eV. This implies that the source spectrum should be 
roughly N{E) cx E'" with a = 2.3 - 2.4, if we as- 
sume an energy- dependent path length, A{R) cx 
(where R = pc/Ze is the rigidity) (Ave et al. 2009). 
If in fact the CR source spectrum at SNRs, G{p), is 
assumed to be released into the ISM at the end of ST 
stage, N{E)dE cx G{p)p^dp is too flat to be consistent 
with the observation. Thus from the spectra Gij)) in 
Fig. 5 we can infer that SNRs expanding into warm or 
intermediate phases of the ISM cannot be the dominant 



sources of Galactic CRs. 

Even with the hot ISM models, the canonical test- 
particle spectrum, N{E) cx E~^ might be still too flat. 
If we consider the effects of Alfven wave drift, however, 
the modified test-particle slope will be given by Eq. 
(2) for strong plane-parallel shocks. One can estimate 
that VA w lOOOkms-i for uh = O.OOScm^^ and Bq = 
30/iG, which leads to a « 2.3. We show in Fig. 6 the 
CR spectra for inefficient injection models with cb = 
0.2. The spectra are less flat, compared with those of 
efficient injection models shown in Fig. 5. Especially, 
the HA model with = 0.2 has the slope, a = 2.1 — 2.3 
for 10^"'^ < E < 10^^ eV. This could be more compatible 
with observed J{E) at Earth. 

Fig. 7 shows the spectra for the models with injec- 
tion recipe B (-Rinj = 3.8). Again G{p) of HB model 
shows the slope, a = 2.1 - 2.3, for 10" <E < 10^^ eV. 
In fact these spectra are quite similar to those for HA 
model shown in Fig. 6. 
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(c) Energy Conversion Factor 

Finally, Fig. 8 shows the integrated energies, Ei/Eo = 
47r J Cir^dr, where eth, ekin, and ecfl are the densities 
of thermal, kinetic and cosmic ray energy, respectively. 
The kinetic energy reduces only slightly and is simi- 
lar for all models. The total CR energy accelerated by 
t/to = 15 is Eck/Eo = 0.35, 0.20, and 0.05 for WA, 
MA, and HA models, respectively, for es = 0.2. In the 
efficient injection models with cb = 0.25 or i?inj = 3.6, 
the evolution of SNRs is quite similar, and the CR en- 
ergy fraction approaches to Ecr/Eo = 0.56, 0.43, and 
0.25 for WA/WB, MA/MB, and HA/HB models, re- 
spectively. So in terms of the energy transfer fraction, 
the CR acceleration in the warm ISM models seems 
to be too efficient. But one has to recall that the CR 
injection rate may depend on the mean magnetic field 
direction relative to the shock surface. In a more realis- 
tic magnetic field geometry, where a uniform ISM field 
is swept by the spherical shock, only 10-20 % of the 



shock surface has a quasi-parallel field geometry (Volk 
et al. 2003). If the injection rate were to be reduced sig- 
nificantly at perpendicular shocks, one may argue that 
the CR energy conversion factor averaged over the en- 
tire shock surface could be several times smaller than 
the factors shown in Fig. 8. 

On the other hand, Giacalone (2005) showed that 
the protons can be injected efficiently even at perpen- 
dicular shocks in fully turbulent fields due to field line 
meandering. In such case the injection rate at per- 
pendicular shocks may not be much smaller than that 
at parallel shocks and the CR energy conversion may 
be similar. Then SNRs in the warm phase of the ISM 
seem to generate too much CR energy. In order to meet 
the requirement of 10 % energy conversion and at the 
same time to reconcile with the CR spectrum observed 
at Earth, SNRs expanding into the hot phase of the 
ISM should be the dominant accelerators of Galactic 
CRs below lO^^eV. 
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V. SUMMARY 

The evolution of cosmic ray modified shocks depends 
on complex interactions between the particles, waves 
in the magnetic field, and underlying plasma fiow. We 
have developed numerical tools that can emulate some 
of those interactions and incorporated them into a ki- 
netic numerical scheme for DSA, CRASH code (Kang et 
al. 2002, Kang & Jones 2006). Specifically, we assume 
that a Bohm-like diffusion arises due to resonant scat- 
tering by Alfven waves self-excited by the CR stream- 
ing instability, and adopt simple models for the drift 
and dissipation of Alfven waves in the precursor (Jones 
1993; Kang & Jones 2006). 

In the present paper, using the spherical CRASH 
code, we have calculated the CR spectrum accelerated 
at SNRs from Type la supernova expanding into a uni- 
form interstellar medium. We considered different tem- 
perature phases of the ISM, since the shock Mach num- 
ber is the primary parameter that determines the ac- 
celeration efficiency of DSA. One of the secondary pa- 
rameters is the fraction of particles injected into the 
CR population, ^, at the gas subshock. Since detailed 
physical processes that governs the injection are not 
known well, we considered two injection recipes that 
are often adopted by previous authors. 

The main difference between the two recipes is 
whether the ratio of injection momentum to thermal 
peak momentum, i.e., Pinj/pth, is constant or depends 
on the subshock Mach number. It turns out the CR 
acceleration and the evolution of SNRs are insensitive 
to such difference as long as the injection fraction is 
similar. For example, the models with injection recipe 
A with €b = 0.23 and the models with injection recipe 
B with i?inj =3.6 show almost the same results with 
similar injection fractions, ^ « 10~^'^ — 10"'^. 

In general the DSA is very efficient for strong SNR 
shocks, if the injection fraction, ^ ^ 10~^'^. The CR 
spectrum at the subshock shows a strong concavity, 
not only because the shock structure is modified non- 
linearly by the dominant CR pressure, but also be- 
cause the SNR shock slows down in time during the 
ST stage. Thus the concavity of the CR spectrum in 
SNRs is more pronounced than that in plane-parallel 
shocks. However, the volume integrated spectrum, 
G{p), {i.e., the spectrum of CRs confined by the shock 
including the particles in the upstream region) is much 
less concave, which is consistent with previous stud- 
ies {e.g., Berezhko & Volk 1997). We have shown 
also that G{p) approaches roughly to time-asymptotic 
states, since the CR pressure decreases as t~^^^ while 
the volume increases as R^r^ cx t^/^. This in turn makes 
the total CR energy converted {Ecr) asymptotes to a 
constant value. If we assume that CRs are released at 
the break-up of SNRs, then the source spectrum can be 
modeled as N{E)dE = G{p)p^dp. However, it is a com- 
plex unknown problem how to relate G{p) to the source 
spectrum N{E) and further to the observed spectrum 



J{E). 

In the warm ISM models (Tq = 3 x lO^'K, nu = 
0.3cm~^), the CR acceleration at SNRs may be too ef- 
ficient. More than 40% of the explosion energy {Eo) 
is tranferred to CRs and the source CR spectrum, 
N{E) cx E~°' with a k, 1.5, is too flat to be consis- 
tent with the observed CR spectrum at Earth (Ave et 
al. 2009). In these models with efficient injection and 
acceleration, the flow structure is significantly modified 
with P2/P0 ~ 7.2-7.5 for WA/WB models. 

In the intermediate temperature ISM models {Tq = 
lO^K, njj = 0.03cm~'^), the flow structure is still signif- 
icantly modified with P2/P0 ~ 5.7-6.0 and the fraction 
of energy conversion, Eck/Eq k. 0.2 - 0.4 for MA/MB 
models. 

Only in the hot ISM model (Tq = lO'^K, Ur = 

0.003cm~^) with inefficient injection (e^ = 0.2 or 
-Rinj > 3.8), the shock structure is almost test-particle 
like with P2/P0 ~ 4.2-4.4 and the fraction of energy con- 
version, Ecr/Eq « 0.1 - 0.2 for HA/HB models. The 
predicted source spectrum G{p) has a slope q = 4.1—4.3 
for 10-^-^ < E < 10^^ eV. Here drift of Alfven waves rela- 
tive to the bulk flow upstream of the subshock plays an 
important role, since the modified test-particle slope, 
Qtp = 3(uo — va)/ {uq — VA — U2), can be steeper than the 
canonical value oi q = 4 for strong unmodified shocks. 
With magnetic fields of Bq = 30pG, the Alfven speed 
is VA ~ lOOOkms"^, and so the modified test-particle 
slope is a ~ 2.3. This may imply that SN exploding 
into the hot ISM are the dominant sources of Galac- 
tic CRs below lO^^eV. One might ask if the magnetic 
field amplification would take place in the case of such 
inefficient acceleration, since the magnetic field energy 
density is expected to be proportional to the CR pres- 
sure. An alternative way to enhance the downstream 
magnetic field was suggested by Giacalone & Jokipii 
(2007). They showed that the density fluctuations pre- 
existing upstream can warp the shock front and vor- 
tices are generated behind the curved shock surface. 
Then vortices are cascade into turbulence which am- 
plifies magnetic fields via turbulence dynamo. 

Finally, in all models considered in this study, for 
Bohm-like diffusion with the amplified magnetic field in 
the precursor, indicated by X-ray observations of young 
SNRs, the particles could be accelerated to iJmax ~ 
10^^-^ZeV. The drift and dissipation of faster Alfven 
waves in the precursor, on the other hand, soften the 
CR spectrum and reduce the CR acceleration efficiency. 
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